Coulomb drag in graphene near the Dirac point 
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We study Coulomb drag in double-layer graphene near the Dirac point. A particular emphasis 
is put on the case of clean graphene, with transport properties dominated by the electron-electron 
interaction. Using the quantum kinetic equation framework, we show that the drag becomes T- 
independent in the clean limit, Tt — > oo, where T is temperature and 1/r impurity scattering rate. 
For stronger disorder (or lower temperature), Tt <C 1/a 2 , where a is the interaction strength, the 
kinetic equation agrees with the leading-order (a 2 ) perturbative result. At still lower temperatures, 
Tt <C 1 (diffusive regime) this contribution gets suppressed, while the next-order (a 3 ) contribution 
becomes important; it yields a peak centered at the Dirac point with a magnitude that grows with 
lowering Tt. 

PACS numbers: 72.80.Vp, 73.23.Ad, 73.63.Bd 



Frictional drag in double-layer systems consisting of 
two closely spaced, but electronically isolated conductors 
is a well established experimental tool for studying the 
microscopic structure of solids [3-13]. In such an experi- 
ment a current I\ is passed through one of the conductors 
(the "active" layer) and the induced voltage drop V2 is 
measured along the other ("passive") layer. The ratio of 
this voltage to the driving current po = —V2/I1 (known 
as the drag coefficient or the transresistivity) is a measure 
of both the inter-layer interaction [l|, Q and the micro- 
scopic state of the layers. At low temperatures the 
drag effect is dominated by direct Coulomb interaction 
between the carriers in the two layers. 

The physics of Coulomb drag is well understood if both 
layers are in the Fermi liquid state 0, EH- The Me- 
tric field in the passive layer is induced by exciting pairs 
of electron-like and hole-like excitations in a state with 
finite total momentum. The momentum is transferred 
from the current-carrying state in the active layer by the 
inter-layer Coulomb interaction. The inter-layer momen- 
tum transfer can be described by the effective relaxation 
rate t7} . The most basic qualitative features of the drag 



measurement P, [l(J II | can already be inferred by esti 
mating r^ 1 with the help of Fermi's golden rule, where 
it is crucial to take into account the energy dependence 
of the density of states (DoS) and/or diffusion coefficient 
D: indeed, the current-carrying states can be charac- 
terized by non-zero total momentum only in the case of 
electron-hole asymmetry. 

The drag coefficient pu and momentum relaxation rate 
t^, 1 can be related using a simple Drude-like model. Con- 
sider the phenomenological equations of motion, assum- 
ing for simplicity that both layers are characterized by 




FIG. 1: (Color online) Drag coefficient in the ballistic regime 
as a function of carrier densities (in units of 10 11 cm" 2 ) for 
d = 9 nm. The left panels show po at T — 250 K, with the 
upper panel corresponding to ultra-clean graphene r _1 = 0.5 
K and the lower left panel showing the evolution of po with 
increasing disorder from r" 1 = to r _1 = 50 K. The right 
panels show po for r _1 = 50 K. The four curves on the lower 
panel correspond to T = 150 K, 200 K, 250 K, and 300 K. 
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where J 1(2) is the average current density in the active 
(passive) layer, E\(2) is the electric field in the two layers, 
and r is the impurity scattering time. Noting that in the 
drag measurement no net current is allowed to flow in the 
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passive layer j 2 = 0, we arrive at the Drude-like formula 
Pd = -p\i = (e 2 nT D /m) 1 . (2) 

Combining Eq. (J2J with the Fermi's golden rule estimate 
for Tjj 1 one can estimate the drag coefficient. More rig- 
orous calculations based on either the diagrammatic per- 
turbation theory [10( or the kinetic equation [lffl confirm 
the "Fermi-liquid" result 
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where Mi (2) is the chemical potential of the active (pas- 
sive) layer and Ayi is determined by the matrix elements 
of the inter-layer interaction (the precise form of A 12 as a 
function of the inter-layer spacing d depends on whether 
transport in the two layers is ballistic or diffusive jiol]). 

Even though the drag coefficient (|3|) is apparently in- 
dependent of the impurity scattering time r, tr ansp ort 
properties of each individual layer are usually [3, [l(| as- 
sumed to be dominated by disorder, r <C td- In particu- 
lar, solving Eq. ([1} for the resistivity one finds the usual 
Drude formula. In contrast, the behavior of clean double- 
layer systems, i.e. with r ^ td, is less trivial. In this 
case, the last term in Eq. (JXJ) may be neglected leading 
to the non-zero result for the single-layer resistivity 



pil = -pi2 = (e 2 nT D /m) = p D . 



(4) 



Note, that the system is still characterized by the infinite 
conductivity (p _1 = 00), as expected for disorder- free 
conductors on the grounds of Galilean invariance. 

The physical picture of the drag effect outlined so far 
is based on the following assumptions: (i) each of the 
layers is assumed to be in a Fermi-liquid state, which at 
the very least means Pi(2) 3> T; (ii) electron-electron in- 
teraction does not contribute to the transport scattering 
time; (iii) the inter-layer Coulomb interaction is assumed 
to be weak enough, a — e 2 /(Hvf) <C 1, such that pp is 
determined by the lowest-order perturbation theory [10j . 

Lifting one or more of the above assumptions leads to 
significant changes in the drag effect 0-Q ■ In this 
Letter we focus on the system of two parallel graphene 
sheets 0-S El-Ell , which offers a great degree of control 
over the microscopic structure of the two layers. Indeed, 
using hexagonal boron nitride as a substrate 0, H2| , one 
can decrease disorder strength in the system and reach 
the regime, where transport properties of the two layers 
are dominated by electron-electron interaction, r 3> r ee . 
Moreover, the carrier density can be electrostatically con- 
trolled allowing one to scan a wide range of chemical po- 
tentials from the Fermi liquid regime to the Dirac point. 

While inapplicable to massless fermions in graphene, 
the equations of motion ((I) provide an expectation of 
non-zero resistance in the case of the ultra-clean system. 
Below, we use the quantum kinetic equation (QKE) ap- 
proach 23|, [24 1 to derive hydrodynamic equations |25[ 
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FIG. 2: (Color online) Drag coefficient in the case of identical 
layers in different parameter regimes for p <C min(T/a,v/d). 
The bottom row of results (for r _1 <C a 2 T and below the 
curve 2, t' 1 < a 2 T 2 /p) are obtained from the solution 
of the QKE (O. The curve 1 (r _1 = a 2 p 2 /T) separates 
the two regimes in Eq. (|15[) . The middle row of results (for 
a 2 T r -1 <C T) corresponds to the region where the results 
and the applicability of the QKE overlap with those of the 
perturbation theory of Ref. [l6| (for p 2> T the results above 
and below the curve 2 contain different numerical factors). 
The third-order contribution p^y = 0(a 3 ) resulting in small 
non-zero drag at p = is shown in red. The upper row of 
results (t -1 3> T) corresponds to the diffusive regime [see 
Eqs. (Tj7|) and |[T8]>]. where p^ saturates for r _1 > T/a 2 ). 



that generalize Eq. ([T]) for interacting Dirac fermions 
in graphene. Solving these equations (or equivalently, 
the QKE) we confirm that the system of two ultra-clean 
graphene sheets is indeed characterized by a non-zero, 
but degenerate resistance matrix whose elements satisfy 
Eq. (HJ), with pa shown in Fig. Q] 

Kinetic equation. — We now briefly outline the deriva- 
tion of the QKE for double-layer graphene structures and 
its solution in the ballistic regime (see Supplemental Ma- 
terial 27]). Consider an infinite sample in an infinitesi- 
mal, homogeneous electric field E\ applied to the active 
layer. The response of the system to the field can be de- 
scribed by the small non-equilibrium corrections /ii(2) to 
the Fermi distribution functions defined by 



n<(e,v)=nj?(e)+T 



de 
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where the eigenstates of the Dirac Hamiltonian H = vcrp 
are labeled [271 ] by their energy e and the velocity unit 
vector v; the momentum of the particle is p = ev/v. 

Small corrections hu2) can be found by linearizing the 
QKE [26| 
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where the linearized pair-collision integrals are given by 
hj = - J d2d3d4 II".://,. i - h i>2 + hj, 3 - h jA ), 

W tj = S(pi - p 2 +P3 - Pi) S(ei - e 2 + £3 - (a) 
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vii \ttH( ,,2! + viv 2 1 +v 3 v 4 

#1,2:3,4 = \U 3 {PI -pa)\ = = , 



two-particle scattering 1 — > 2 and 3^4 and the corre- 
sponding Dirac factors. Here we take into account only 
the Hartree interaction term: there is no exchange inter- 
action between the layers, whereas within the layers the 
Hartree term dominates in the large- N limit (N is the 
number of electron flavors; physically, N — 4 due to spin 
and valley degeneracy). 

The peculiarity of the inelastic scattering in the Dirac 
spectrum is two- fold. First, since the velocity v = v 2 p/e 
is independent of the absolute value of the momentum, 
total momentum conservation does not prevent velocity 
(or current) relaxation. As a result, the intralayer col- 
lision integral Zy yields a non-zero transport relaxation 
rate due to electron-electron scattering. 

Second, the scattering of particles with almost collinear 
momenta is enhanced since the momentum and energy 
conservation laws coincide for collinear scattering. This 
restricts the kinematics [23|, 0, [28| of the Dirac fermions 
leading to the singularity in the collision integral. This 
singularity leads to the fast thermalization of particles 
within a given direction, which justifies the Ansatz: 



hi(e,v) 



X«+X« e/T)eEv/T\ 



(9) 



The Ansatz ^ retains the only two modes for which 
the collision integral is not singular: the "momentum 

mode" Xv\ which nullifies the collision integral due to 
momentum conservation, and the "velocity mode" Xv\ 
which nullifies 7^ in the case of collinear scattering. The 
same kinematic restrictions lead to fast uni-directional 
thermalization between the layers. This allows us to set 



(2) - 

, and hence reduce the QKE for the double- 



layer setup to a 3 x 3 matrix equation. 

Consider for simplicity the case of identical layers (for 
the more general case of \i\ 7^ /12 see Supplemental Mate- 
rial 27|). Integrating the reduced QKE over the energies, 
we arrive at the set of steady-state hydrodynamic equa- 
tions in terms of the particle currents 



J, = -NT 
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and the total momentum P = etQC\(E\ + E^t: 
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where l ee (D) = [(<5o + °i)Cl + 2a ^x)C2]/T ee ( D ), the 
intra- and inter-layer electron-electron transport scatter- 

2|V ^ //r,rp\l\ 
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and Cfc are the Pauli matrices in "layer space" . The co- 
efficients Ci(2) represent the average energy and energy 
variation, while £q = 2TJ{1}/N is a typical energy: 
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The hydrodynamic equations (1111) generalize the equa- 
tions of motion (TT|) to the case of Dirac fermions in 
graphene. The kinematic peculiarity of Dirac fermions 
manifests itself in the appearance of the total momen- 
tum, which entangles the electric fields in the two layers. 
Solving the hydrodynamic equations (|llj) we find 
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KC 2 {TT D )^+Cl[T- e 
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(13) 



For a clean system, the resistivity matrix is degenerate 
and the drag coefficient is given by 



Pd(t 



= (h/e 2 )(C 2 /e Q ) 



D 



(14) 



which remains non-zero po ~ (H/e 2 )a 2 even at the Dirac 
point fi — 0, where it is determined by t~} ~ a 2 T (to the 
second order in the inter-layer interaction T I ^ 1 (fi = 0) = 
0, while the third-order contribution r z ^ 1 (/i = 0) ~ a 3 T 
is subleading; the latter is expected to dominate the effect 
for sufficiently strong disorder, see below). 

Equation (Ti"3"|) gives the general expression for the drag 
coefficient in the ballistic regime based on the solution of 
the QKE ©. For arbitrary parameter values this ex- 
pression is to be evaluated numerically (see Fig. Q] for 
the numerical results) . Analytical expressions can be ob- 
tained for various limiting cases (summarized in Fig. [2]) . 
Below, we discuss the asymptotic behavior of pr> for the 
case of two inequivalent layers [13] focusing on the ex- 
perimentally relevant case [7H9( Td/v < 1 and analyzing 
the evolution of po with increasing disorder strength. 
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Ballistic regime. — For weak disorder a 2 Tr ^> 1 (or 
r -1 <C r^ 1 ) and neglecting the third-order contribution 
to Tjy 1 , we find for pr> near the Dirac point 



Pd{p% < T) 



2.87^ cf-x- 

e 2 /if 



^1 



0.49T/(a 2 r)' 



(15) 



where r^ 1 ~ a 2 P1P2/T , t^} ~ a 2 T, eo ~ T, and C2 ~ 1. 
The value of p£> precisely at the Dirac point depends 
on the experimental set-up. For clean samples, if one of 
the chemical potentials remains non-zero, while the other 
is scanned through the Dirac point 0, then pd(pi = 
0, p% =7= 0) = 0, similar to Ref. H|. On the contrary, if 
both chemical potentials are driven through the Dirac 
point simultaneously then Eq. (|15[) predicts a non- 
vanishing value of pd{p-i — ±P2 = 0) 7^ 0, see Fig. [TJ 

For intermediate disorder strength a 2 T <C t _1 <C T 
the applicability region of the QKE overlaps with that 
of the conventional perturbation theory developed in 
Ref. and we recover perturbative results, see Fig. 

For even stronger disorder (or at low temperatures) 
Tt <C 1 the electron motion becomes diffusive. In this 
case the kinematic restrictions are relaxed and the Ansatz 
(J9j> is no longer justified. However, in this regime, the 
perturbative approach is applicable and allows for a stan- 
dard description of the diffusive transport. 




FIG. 3: (Color online) Schematic view of the drag coefficient 
at low temperatures: the second-order contribution (solid 

line) and the third-order contribution pf/ (blue dashed line). 
The arrows indicate the tendency of the two terms with the 
decrease of temperature T — > 0. 

Diffusive regime. — The lowest-order perturbative cal- 
culation [lOj amounts to evaluation of the Aslamasov- 
Larkin-type diagram for the drag conductivity given by 
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where "Dyi is the retarded propagator of the inter-layer in- 
teraction and r"(w, q) is the non-linear susceptibility [in 
fact, all previous studies of the Coulomb drag in graphene 
12j-|20( focused on Eq. fllS])]. In the diffusive regime, 



r™(u;, q) can be found using the Ohm's law and the con- 
tinuity equation [3(| T = eq(pa j dn)hxJl R . All micro- 
scopic details are now encoded in the diffusion coefficient 
and the density dependence of the single- layer conductiv- 
ity a. Close to the Dirac point p <C T <C t _1 the deriva- 
tive da/dn ~ nv 2 T 2 (independently of the precise nature 
of impurities). After this the evaluation of Eq. (fT6|) is 
rather standard (except that, in contrast to Ref. |l0|, the 
Thomas- Fermi screening length is much longer than the 
inter-layer spacing xd Cl) and yields 



pg } {fa « T « r- 1 ) ~ (hle 2 )a 2 p Y p 2 TT Z 



(17) 



This result vanishes at the Dirac point as a consequence 
of the electron-hole symmetry. 

The importance of the electron-hole asymmetry for the 
Coulomb drag follows from Eq. (fT6]): the non-linear sus- 
ceptibility can be thought of as a measure of the asym- 
metry. However, Eq. (|16p is only the lowest-order contri- 
bution to <td- Under standard assumptions of the Fermi- 
liquid behavior in the two layers (p ^> v/d 3> T, pr 1), 
this contribution indeed dominates the observable effect. 
On the contrary, in the vicinity of the Dirac point in 
graphene, the next-order contribution p^ 1 [29[ becomes 
important since it is insensitive to the electron-hole sym- 
metry and thus does not vanish at the Dirac point. 

The explicit results of Ref. [2^ were obtained in the 
usual limit xd 3> 1. Extending these calculations to the 
opposite case xd < 1 we find close to the Dirac point 



(3) 

Pd 



Pi < T < T" 



<a~ 2 T 



(h/e 2 )a 3 (TT)- 3 / 2 , 

(18) 

and p^j ~ h/e 2 for t _1 ~s> a~ 2 T. Away from the 
Dirac point this contribution decays as a function of 
the chemical potential p^\pT ^> max[l, a _1 (Tr) 1 / 2 ]) ~ 
(H/e 2 )(pT)~ 3 and rapidly becomes subleading. As a re- 
sult, p^p is only detectable at low T and p, see Fig. |3] 

While estimating p D at the Dirac point, we assume 
the single-layer conductivity a ~ e 2 /h discarding local- 
ization effects. Indeed, experiments on high-quality sam- 
ples show T- independent a down to T = 30 mK 3l| , that 
can be explained by the specific character of disorder in 
graphene [32j. 

Summary. - We have studied Coulomb drag in 
double-layer graphene structures. By using the QKE for- 
malism we have shown that for weak disorder (or high 
T; ballistic regime) po near the Dirac point is given 
by Eq. (fl"5|) . see also Fig. Q] which is consistent with 
Ref. H. For a 2 TV <C 1, the solution of the QKE agrees 
with the perturbative calculation of Ref. 1Q. For even 
stronger disorder (or lower T; diffusive regime) sublead- 
ing third-order contribution dominates the effect in qual- 
itative agreement with the experimentally observed peak 
at the Dirac point at low temperatures Q. A possible 
alternative origin of the low-T peak at the Dirac point is 
a collective state of the double-layer system, either due 
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to strong interaction [33| (not explored here) or corre- 
lated disorder 34|, |35[ (discussed in Supplemental Mate- 
rial [13). 
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Coulomb drag in graphene near the Dirac point: Supplemental Material 
1. Kinetic equation approach 

In this section we provide details of derivation of the resistivity tensor in double layer graphene from the kinetic 
equation approach. 

A. Kinetic equation 

Eigenstates of the massless Dirac Hamiltonian H = vcrp are characterized by the values of momentum p and the 
discrete variable x = il indexing conduction and valence bands. In this representation, energy and velocity are 
e = av\p\ and v = x v p/p- It is, however, more convenient to label the eigenstates by their energy e and the unit 
velocity vector v. The momentum of the particle is then p = e\r/v and the normalization of the states reads 



\e\dedv 

|e,v)(e,v| = 1. 



(2ttw) 2 



(19) 



We consider an infinite double-layer graphene sample (layers 1 and 2) in a homogeneous electric field Ei applied 
to the active layer 1. Assuming weak electric field, we start with the linearized kinetic equation: 

dhi eE lV ht 

-^r + ~tf~ = \-hi{hi} + Ji 2 {/H,M> 

at 1 t 



^ = -— + I 2 2{h 2 } + I 2 i{h 2 ,h 1 }. 
at t 

Here the nonequilibrium correction hi to the Fermi distribution function is defined by 

7ii(e, v) = n F {e) + T — hi{e, v), 



ch 



(20) 



(21) 



and Iij is the linearized pair-collision integral: 



Iij = -N / de 2 de 3 de 4: / d\r 2 dw 3 dv i v(e2)v(e 3 )v(e 4 )W l: > (1, 3; 2, 4) [(hi(e u vi) - hi(e 2 , v 2 ) + hj(e 3 , v 3 ) - ^(e 4) v 4 )] , 



W ij {l, 2; 3, 4) = S( Pl - p 2 + p 3 - p 4 ) 6{ex - e 2 + e 3 - e 4 ) 



cosh 



2T 



2 cosh cosh cosh ^p- 



(22) 

— jr«(l,2;3,4). (23) 



Here v(e) is the density of states for one of N flavors (per spin and per valley in graphene, where N = 4). We assume 
formally the large N limit and neglect the intralayer exchange interaction. We further assume that the scattering does 
not mix flavors (i.e., we neglect the intervalley scattering due to Coulomb interaction): states 1(3) and 2(4) belong 
to the same flavor, which gives the overall factor TV. The kernel 



(24) 
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contains the interaction matrix element My describing the collision of two particles 1 — > 2 and 3 — > 4 and the 
corresponding Dirac factors. Within the Golden-rule approximation, this matrix element is given by the Fourier 
component of the interaction potential: 



where 



with 



Mlp=U§ ) (px-pa), (25) 



U i0) (q) = V (q)[J qd C ^ d ), (26) 



2?re 2 

V Q (q) = ^-. (27) 



Further, one can generalize the collision integral to the case of the RPA-screened interaction. Then 



where 



Mij = l/Jf A ( Pi - p 2 ,v Pl - v P2 ), (28) 

Mq) 

[1 + V (g)IIi (g, w)] [1 + U (q)n 2 fa w)] - e^ d V 2 (g)IIi (q, lo)I1 2 (q, u) 

(29) 



l + V (q)n 2 {q,L})(l~e- 2c i d ) e"« d 

e~« d l + Vb(q)ni(< ? ,w)(l-e- 2 « <i ) 

where w) is the polarization operator in layer i. 

We will focus on the experimentally relevant case of closely located layers, Td/v 1. Furthermore, here we will 
restrict our consideration to the case of relatively low concentrations, such that fid/v <C 1 (the situation with large 
interlayer distance will be considered in detail elsewhere) . Under these conditions we can set d = in the interaction 
matrix elements so that the intralayer and interlayer interactions are just the same. We further assume that the 
interaction coupling constant is small 

e 2 

a = — < 1. (30) 

v 

For simplicity, we treat impurity scattering within the relaxation time approximation with an energy-independent 
transport time r. Generalization to the more realistic case of Coulomb impurities with an energy-dependent transport 
time will be discussed elsewhere. 



B. Collinear-scattering singularity 

The momentum and energy conservation establishes severe kinematic restrictions on the scattering in systems with 
linear spectrum 23[ . This can be easily seen, when one rewrites the product of delta- functions in Eq. ([2"B")) as integrals 
over q = Pi — P2, ui — e\ — e 2 , and two angles (pits) between q and Pi(3), respectively: 

f d 2 q 

d>(viei - v 2 e 2 + v 3 e 3 - v 4 e 4 ) = J £(viei - v 2 e 2 - q)^(v 3 £3 - v 4 e 4 + q), 

6 (ex - e 2 + £3 - £4) = / dujS(e 1 - e 2 - u)5(e 3 - e 4 + w) 3 
J 00 

(31) 

This allows one to integrate out p 2 and p 4 (v 2>4 and e 2i4 ) in the collision integral, leading to the product 



S(ei - oj - ^ (\ + q 2 v 2 - 2eiqvcos4>i)S(e 3 +w — \Je%+ q 2 v 2 + 2e 3 gwcos0 3 ), (32) 
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which sets 



gV - uj 2 + 2e x uj 

ei-w = upi-q =3- cos 0i = , (33) 

zeiqv 

cj 2 - c?V + 2e 3 u 

e 3 +u;=|p 3 + q| =4> cos0 3 = . (34) 

2e 3 gw 



When calculating the scattering rates using the collision integral (|22| . the angular integration over cf>i and 3 removes 
the delta-functions, producing the factor 

e x -uj e 3 +uj 4(ei - w)(e 3 + uj) . . 

(35) 



ei^|sin</)i| e 3 gt;|sin^ 3 | ( q 2 v 2 ~ U j 2 )y / [(2e 1 - uj) 2 - q 2 v 2 ][(2e 3 + uj) 2 - q 2 v 2 ] 
In combination with the Dirac factors from Eq. ([24|). 



this yields 



H _i_ uu ^ (2 £l -^) 2 - g V (2 £3 +^) 2 -g 2 t ; 2 

(1 + viv 2 )(l + v 3 v 4 ) = — r — - — — , (36) 

2e 1 {ex-uj) 2e 3 (e 3 + u;) 



,/(2e 1 -uj) 2 -q 2 v 2 ^(2e 3 +uj) 2 -qh? 



q 2 v 2 — uj 2 ei e 3 



Therefore, further integration over q or uj generically produces a logarithmic divergence at uj — ±gu, which stems 
from the collinear scattering <fi 3 = <\>\ = or tt, see Eq. (|34|) at the light cone. This divergence reflects the fact that 
for a linear spectrum the momentum and energy conservation laws coincide in the "one-dimensional" collinear case. 
Note that this enhancement of the collinear scattering is not restricted to the case of undoped graphene. 

In order to regularize this divergence, one has to go beyond the Golden-rule level and take into account the 
screening of the interaction (which in the clean case is perfect exactly on the light cone) and renormalization of the 
spectrum due to interaction (leading to nonlinear corrections). These mechanisms lead to the appearance of a large 
factor | ln(a)| 3> 1 in generic relaxation rates in graphene. In disordered graphene, this singularity is also cut off by 
disorder-induced broadening of the momentum-conservation delta-function. 



C. Ansatz 



The singularity in the collinear scattering in graphene leads to the fast thermalization of particles within given 
direction within each of the layers. Clearly, in an external electric field E, the linearized nonequilibrium correction to 
the distribution function is proportional to the driving term Ev. Therefore, the nonequilibrium correction 

hi(e,v) = x(t)^r 

is characterized by some function of energy x( e )- O ne can formally expand this function in e/T: 

oo 

X( e ) = $>„(e/2T. 

71 = 

The action of the collision integral (which contains the combination hi — hi + h 3 — h^) on this function generates the 
equilibration rate in all terms except for n — 0,1. Indeed, for n — 0, the combination 

^-^+4°)-/ 1 foc X0 (v 1 -v 2+ v 3 -v 4 ) 

contains the differences of velocities. This cancels the kinematics-induced divergency tx |vi — v 2 + v 3 — For 
n = 1 the combination 

h^ - hip + /i 3 1) - oc xi(eivi - e 2 v 2 + e 3 v 3 - e 4 v 4 ) = Xi(Pi P2 + P3 P4) (38) 

contains the change of the total momentum of two colliding particles, which is exactly the argument of the momentum- 
conservation delta-function. All other contributions with n > 1 produce a relaxation rate enhanced by the collinear 
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scattering. The corresponding values of \ n are therefore strongly suppressed compared to xo and xi, which justifies 
the following Ansatz 231 ] : 



11 \ ( i i e ^f I \ eEv ( i i e \ eEv ( i i e \ 
ft»(e,v) = I X M + Xt jl I -ya" = (^Xo + Xi^ J -yj- = (X„ + X P ^ J 



e \ eEv 
T/ T 2 ' 



(39) 



This correction to the distribution function contains only the two modes (proportional to velocity and momentum and 
characterized for each layer by the two constants Xv = Xo and Xp = Xli respectively), that nullify the collision integral 
in the case of collinear scattering. The notation Xfj. and XT is chosen to emphasize that, after the linearization with 
respect to E, these quantities reflect the angular-dependent corrections to the chemical potential and temperature, 
respectively, in the direction-equilibrated distribution function: 



n(e,v) = 



1 



1 + cxp 



M(v) 



n F (e) - 



2T(v) . 
1 dnp(e 
2T de 



n F (e) 



dn F {e) 
de 



2T 



2T 2 



{[^(v)-^T(v)] +^T(v)} 



(40) 



(41) 



The Ansatz Eq. (|39p greatly simplifies the solution of the kinetic equation, replacing the integral equation by a 
matrix one. Furthermore, the same kinematics-induced singularity in the collinear scattering as in la appears in also 
the intralayer collision integrals 1^. This implies fast unidirectional thermalization between the layers. Therefore, 

we set Xt^ = Xt\ anc ^ hence reduce the kinetic equation for the double-layer setup to a 3 x 3 matrix equation 
("three- mode approximation"). 



D. Hydrodynamic equations 

In each of the two layers, we introduce the particle currents (the total velocities) 

r 5 n W r 

J* = -NT / de f(e)-gf- / dvvh a {e,v), 

and energy currents (or, equivalently, the total momenta) 

f dn (i) f 

Pi = -AM de v(e) e / dv vhi{e,v), 

and substitute the Ansatz, Eq. [3H into these expressions, which yields 



P, 



N 
~2 

N 

~2 



/|(*) V W i 4(') v (0 
A o Xo + A i Xi 

A i Xo + A 2 Xl 



E 
E. 



Here 



v 2 f , , . drip ( e \ » 

t J de ^ "flf (r) ' 



In terms of the currents, the fast interlayer thermalization (Xt = Xt J translates into the relation 



B b 2 (P a - B«J a ) = B a 2 {V b - Bjj 6 ), 



where 



r>a _ \a 



B? = 4^, B% = A%- 

^-0 -^0 



(42) 

(43) 

(44) 
(45) 

(46) 

(47) 
(48) 
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The asymptotics of the functions B$ (/Xj/2T) read: 
B«(z«l) = 



ln2/7r n = 
4a; n = 1 , 
9C(3)/2tt n = 2 



B n (a; > 1) = 



| a?| / vr n = 
2x n = 1 
7r|x|/3 n = 2 



(49) 



It is also convenient to define B 2 — B^ + B^ ■ Finally, we introduce the total momentum (total energy current) 

P = P a + P 6 , (50) 

which is not affected by electron-electron collisions due to total momentum conservation in the e-e collision integral. 

These transformations allow us to rewrite the matrix kinetic equation in the "hydrodynamic form" . Here we can 
also introduce electric fields in both layers without doubling the number of relevant modes. Integrating the reduced 
matrix kinetic equation over the energy with 



?(2) 



and 

-iV 



/^ (() ^,.,. w / il£1 _^L,., 



(51) 



(52) 



yields the following steady-state equations (in order to avoid confusion in indices, from now on we label the layers by 
a and b) 



1 



--B«Bt-)J a + \ B\B\-- 



b\2 



Eh 



TD 



_L -B«Bl— 



1 



{Bl) + M 



— J>J„ 

TD 



2 

N~ 



—TB%eE a + I — -]P 



jjTB b eV b 



Bl B\ 

Tee T D 

B 2 



B a 



a\2 



(B?) 



1 

TD 



Tee T D 



Here we have introduced the following effective transport relaxation rates: 



iV 



8T 2 B 2 

TV 
8T 2 B 2 

N 



T D 4T 2 B 2 



d{ei}d{vi} [( Vl - v 2 + v 3 - v 4 ) 2 W QQ + 2 (v x - v 2 ) 2 W ab 
d{ ei }d{vi} [(VX-V2+V3-V4) 2 W bb + 2(v 1 -v 2 ) 2 W i 
fd{ei}d{vi} (vi-v a )(v4-v a ) W ba . 



Here l/r e a e (6) are the intralayer transport relaxation rates describing the velocity relaxation within a layer due to 
inelastic scattering with electrons in the same layer (described by W aa ) and in the other layer (W a& term). The rate 
1/td describes the inter layer velocity relaxation (velocity transfer from one layer to the other due to W ab ): we call it 
the drag rate. The kernels here are related to the kernel of the collision integral ([23]) as follows: 



(53) 

(54) 
(55) 

(56) 
(57) 
(58) 



W y (l,2;3,4) 



cosh 



-W™(1,2;3,4). 



(59) 



IT 
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2. Ballistic regime 
A. Resistivity matrix 

The hydrodynamic equations (|55j) yield the following explicit expressions for the intralayer and interlayer resistivi- 
ties: 



2B 2 



e 2 N{B*) 2 T 



1 



B a 



1 



1 



D 



1 
T 



{Btf 



2B C {B\ 

TD 



IB, 



Pab 



1 

TTD 



+ Bf B\ 



e 2 NBftBfiT i {B af ( B bf 2B a Bl , 



(60) 



(61) 



^~ ^ee Tse 

Pbb = Paa(fl <-> b), and pba = Pab- The drag coefficient is defined as 

PD = -Pab- 

It can also be rewritten in the following form, 



TD 



(62) 



2Bo 



PD 



e 2 NB^B^Ttd 





3) 


(B\ 




\ ' ee 


TD J 


We " 


TD J 




2B\B\ 



TD 



(63) 



which for the disorder-dominated case yields directly the conventional perturbative drag: 

h 2B 2 



For equal layers, we denote 



PD e 2 NB$B%Tt d ' 
e = 2B T/N, C l = B 1 , C 2 =B 2 /B . 



(64) 



(65) 



Then the resistivity matrix reads 

fr_ C2 fl 

e 2 e It 



P 



1 
1 



c 2 



1/r + 2C 2 (l/r ee - l/r D ) 

C 2 C\ (l/r 2 e - 1/r 2 ) 
1/t + 2C? (1/t„ - 1/td) 



l/r ee -1/td 
-1/td l/r ee 




(66) 



Here the first term is the intralayer resistivity determined by disorder, the second term describes the conventional 
Coulomb drag in combination with the intralayer inelastic transport relaxation, and the last term arises due to the 
fast unidirectional thermalization in graphene. 

In the clean limit r = oo the resistivity matrix has the form 



1 



h 2B 2 



J_ / [B\) 2 BIB\ \ 

0/ 



e 2 NT {B af 



2BIB\ 

TD 



r B a\2 



f>a f>b 

B%B\ {Bf) 2 



(67) 



which for equal layers simplifies to 



„ h c 2 ( i i 
p = — — — + — 

e 2 £ \ r ee TD 



1 -1 
-1 1 



(68) 



The off-diagonal component of this matrix is given in Eq. (14) of the main text. 
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B. Asymptotics of the drag coefficient 



The general condition separating the disorder-dominated and Coulomb-dominated transport regimes can be found 
from Eq. (|61[) . where one should compare the two terms in the numerator: 

~ - BtB\r D ( * - -L) . (69) 

' \ 'ee 'ee 1 D / 

In the vicinity of the Dirac point (p a ,b <C T), the intralayer transport relaxation rates are 

~c?NT, (70) 



whereas the drag rate was found in Ref. [l( 



-~a 2 A^, (71) 

T~D 1 



so that 1/Vee S> 1/td. Substituting these results into Eq. (|69|) . we find 

a 2 AT, (72) 



1 flafJ-b T D , 



r 

i.e. the crossover occurs at 1/r ~ l/r ee . 

Away from the Dirac point, T <C A*a,& <C T/a, the drag rate is (for simplicity we set \x a ~ A*b ~ A 4 ) 

1 T 2 u 

a 2 iV— ln^, (73) 

and 

1 1 IT 2 1 
2 « — ■ (74) 

7"ee To ^ee A 4 'ee 

It is worth mentioning that, for /i>T, the intralayer scattering is much less efficient for the relaxation of velocity than 
the interlayer scattering. Indeed, the intralayer (oc W") contribution to l/r ee in Eq. (|58p contains the combination 
of velocities vi — V2 + V3 — V4, which for /1 ^> T is very close to Pi — P2 + P3 — P4- Therefore, at high chemical 
potentials l/r ee is dominated by the interlayer contribution (oc W ah ). Thus the crossover between the two regimes 
occurs for /x>Tat 

T VT ee To/ A* T T ee 

The drag coefficient in the disorder-dominated regime coincides with the perturbative result, 

h 2B 2 



Pd = 



In the opposite, ultraclean case, assuming for simplicity equal layers, we get from Eq. ([68 

h C 2 ( 1 1 \ h 2C 2 



e 2 NB^B^Ttd 

h 2 T 2 (n a +n b ) min{A<a,Ai&} ™ „ „ ^ 1 2)xr T2 

VHaPb) 1 (XT fl 



P D - 2 \ • 1 — 2 

fr T 2 u T 1 T 2 

~4a 2 ±-ln^, T «,;«-, i«a 2 7V— . (77) 

The difference between the disordered (perturbative) and ultraclean (equilibrated) results for the drag is in the 
presence of l/r ee in the latter case. In particular, for equal layers, this enhances the drag by a factor of 2. Since 
these results are qualitativel y th e same, we do not analyze the behavior of the drag at yet higher chemical potentials, 
referring the reader to Ref. [16J. 
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C. Third-order drag rate 

It is important that the second-order (Golden-rule) drag rate vanishes at the Dirac point due to the particle-hole 
symmetry. However, the particle-hole symmetry does not affect the odd-order drag rates [2y]- Near the Dirac point, 
such rates should not depend on /i and hence are proportional to T. Taking into account the second-order matrix 
element ex a 2 , 



M {1 h ] + M { 1 ] 

ab ' ab 



M 



(i) 



ab 



2Re(M« 



M 



(2) 
ab 



(78) 



we estimate 



— ~ a 2 N^—- + a 3 NT, 
t d T 2 



(79) 



where we skip the numerical prefactors in both terms. A similar correction arises in l/r ee , but there it is always 
subleading for a <c 1. Substituting this correction into Eq. (1661) . we find 



PD 





+ ^N 2 


aiT 2 -(a 2 ^+a^ 2 


e 2 NT 1 u 2 N 

+ 

T T 2 


a 2 T+ (a 2 ^+a 3 T^ 







(80) 



(81) 



Clearly, for /i 3> a x ^ 2 T we can disregard the a 3 -terms. In the opposite limit, we neglect the conventional drag term 

h a 3 T + a^ 2 rN 1/2 1 

e 1 1 + a z fi 2 rN t 

Exactly at the Dirac point this yields 



Pd ~ -zct 
e z 



For finite chemical potential, this result is valid for 



T 1 



(82) 



(83) 



In the opposite limit, one can neglect the a 3 -term, yielding 



PD 



h a 4 fi 2 rN h a 2 , 



1/r <C a 2 Nn 2 /T 



e 2 T + a 2 /i 2 TN e 2 a^N^r/T, a 2 Nfi 2 /T < 1/r < aN^/T 



(84) 



3. Diffusive regime 
A. Third-order contribution to the drag 

In this section we analyze the third-order drag [29j j in the diffusive regime Tr <C 1 for the case of high dimensionless 
conductances (per spin and per valley), g = vD ~ > 1, where D = v 2 t/2 is the diffusion coefficient. In this 
limit, one can calculate the prefactor analytically. In the vicinity of the Dirac point /jt < 1, the conductance is 
of order unity. Indeed, experiments on high-quality samples show T-independent a down to T = 30 mK [311 ]. that 
can be explained by the specific character of disorder in graphene |32j |. Therefore, there we also assume the diffusive 
dynamics described by the diffusion propagators 

T>i(q,u>) = - 1 , , qv,uj<^l/T, (85) 
Vi Viq^ — leu 



13 



and polarization operators 



Diq 2 — iuj 



(86) 



where Vi is the density of states of the layer (per spin and valley) i = 1, 2 and Di is its diffusion coefficient. 

The third-order drag resistivity was calculated in Ref. (2l| for the case of large interlayer distance, nd 3> 1. Here 
we generalize this result to the opposite case «d < 1, which is experimentally relevant for graphene near the Dirac 
point. 

The analytic expression for the third-order drag resistivity is given by [29j 

oo 


(2^ J (^ Im [^i(9^)^(<z^)V 12 ( g , w )Vi 2 (|_Q,|_n) Via(| + Q,| + n 



(87) 



Here the thermal factors J-i(w, f2) are given by 



Ji(w, 0) - T— [B(Q + w/2) - - oj/2)] 
oil 



(88a) 



J" 2 (w, Q) = 2 - B(fi + w/2) - B(n - w/2) + B(w), 



(88b) 



where 



BM = ^coth(^) 



3c) 



The propagators of longitudinal vector potentials Vi 2 (q, w) in Eq. (|87|) include the dressing of the vertices by diffusons 



Vi 2 (g, w) 



{Diq 2 — iw)(Z? 2 <7 2 — iw) ' 



(89) 



where the retarded RPA-screened interlayer interaction U^' A (q, uj) is defined in Eq. (l29l) . 

For simplicity, we consider equal layers. For small interlayer distance d <^ vt we have qd <C 1. ft is convenient to 
introduce the inverse screening length 

K = 2ire 2 is, (90) 
where is the thermodynamic density of states per one flavor of particles. Expanding exp(— qd) ~ 1 — qd we get 

(l + U^iqMq^)) 2 - (u$\qMq )U j) 
Nn Dq 2 \ 2 {NKe-i d Dq 2 



1 ne- qd 
v q 

1 K 



q Dq 2 — iuj 
(Dq 2 - iuS) 2 



Dq 2 — iuj 



q [Dq{q + 2Nn) - iuj][Dq 2 {l + Nnd) 



and hence 



Kq 



v [Dq(q + 2Nk) - ioj][Dq 2 (l + Nad) - iuj] ' 
We first consider the case of large interlayer separation, Nnd ^> 1. For Nn ^> max{l/<i, (T/Z)) 1 / 2 } we find 



2Ng Dq 2 Nnd-iuj 



(91) 



(92) 



(93) 



14 



which reproduces the result of Ref. 29] : 



(3) 



h 1 



Pd 



e 2 ^3^3 ( NK d) 2 



For l/d<iVK< (T/D) 1 / 2 , 



and we find 



y> t \ 1 Kq 



l 



v Dq 2 — iuj Dq 2 Nnd — iuj ' 



HI 1 



Pd ~ 72 „3 



e 2 5 3 (Nnd) 2 \ T 

In the opposite case Nnd <C 1 (which is relevant to our problem) , 

v , s 1 ^7 

12 W ' Wj ~ i/ [£>g(g + 2JV K ) - iw] [ZV - iw] ' 

Substituting this into Eq. (j57J, we get 



(94) 



(95) 



(96) 



(97) 



OO 



(3) _ h on^S 2 



(I 



(2tt) 2 J (2tt) 2 
«|q/2-Q| 



1 



Dq 2 — iuj 



Kq 



[Dq{q + 2Nn) - iuj][Dq 2 - iu)] 



[D(q/2 - Q) 2 + 2L>|q/2 - Q,\Nn - i(w/2 - Q)} [D(q/2 - Q) 2 - i{w/2 - O)] 

«|q/2-Q| 

[D(q/2 + Q) 2 + 2£>|q/2 + Q|7Vk - i(w/2 + CI)] [D(q/2 + Q) 2 - i(w/2 + O)] 

The frequency integrals are dominated by ~ 57 ~ T, whereas the momentum integrals are dominated by q 
qr = y/T/D. Therefore, the drag conductivity in Eq. (|98|) can be estimated as 



(98) 



Pd 



(3) 1 1 r j ,.| 



duidQ 



d 2 qd 2 Q 



(Dq 2 + Tf (Dq 2 + T + Dq T N K ) 3 



hi k 

e 2 g 3 \Nk+ y/T/D 



Therefore, for Nk > y/T/D we get 



.(3) 



h 1 



Pd ~ -2 



whilc for TVk < ^T/Z? we find 



e 2 



(3) H 1 fDn 2 ^' 2 
Pd ~ ■> „3 



e 2 g 3 \ T 

In the first case the expression for the third-order drag coefficient is given by 



(99) 



(100) 



(101) 



Pd = | 32T^ / d^*( W> n)Ji(a,,n) / ^ / ^ Im 
o 



1 



Z?g 2 — 
1 



2iVD / ZV - D(q/2 - Q) 2 - j(w/2 - Q) D(q/2 + Q) 2 - i(u/2 + Q) 



(102) 
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The integrals here are now dimensionless [one measures momenta in units of (T/D) 1 / 2 and frequencies in units of T]. 
In the second case the prefactor is again determined by a dimensionless integral: 



o 



1 



Dq 2 — iuj 

q/2-Q| |q/2 + Q| 



[Dq 2 ~iu} 2 [D(q/2-Q) 2 -i{uj/2-n)} 2 [D(q/2 + Q) 2 - i(u>/2 + O)] 2 
For TV <C 1 and /UT 3> 1, we have 



(103) 



k ~ a/x/u, (104) 

which implies 

/I 7 1 a 2 N 2 fi 2 , . 

^» = Vd ° 7 = ^- (105) 

For /iT <C 1, we have 

/c ~ a/vr, (106) 

and hence 
Thus, for 

max{T,a 2 AV/T} «1/t< T/a 2 N 2 

the third-order drag resistivity reads: 

At Tr ~ 1 this result matches the ballistic third-order drag resistivity, Eq. (|82p. 
For yet lower T <C a 2 N 2 /r and /xt <C 1 the third-order drag saturates at 

Finally, for /xt 3> max{l, (Tt) 1 / 2 /^^}, the third-order drag behaves as 

pg) ^ 1 (110) 

4. Correlated disorder 

In the original version of the paper we mentioned the correlations between the disorder potentials of the two layers 



34J as an alternative mechanism leading to a low-T peak in po at the Dirac point. After the submission of the 
original version, we became aware of a preprint by Song and Levitov [35| that focused on such a mechanism. In this 
section we analyze the role of interlayer correlations of disorder potentials (both of short-range and long-range nature) 
within our general framework. This allows us to compare the effect of correlated disorder with the third-order drag 
considered in the main text and in Sections 2C and 3 of the Supplemental Material. 

As emphasized in Ref. [35j |. the correlations between the disorder potentials of the two layers might be especially 
important in drag experiments on graphene near the Dirac point for the two reasons: (i) similarly to the third-order 



drag, it does not require 34| the particle-hole symmetry and hence provides finite drag at the charge neutrality point 
[35j ]: (ii) in contrast to experiments on conventional semiconducting double wells, the interlayer distance in graphene 
experiments is rather small, which enhances the disorder correlations between the layers. In what follows, we analyze 



the two models of correlations: (A) Correlated scattering off common short-range impurities 34| and (B) correlations 
of large-scale inhomogeneities of the chemical potentials in the layers (3|| . 
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A. Short-range correlations: correlated impurity scattering 



Following Ref. [34| . we introduce the matrix of disorder correlators 



(u y 'u 



The values of w^' at 



i ^ j differ from zero due to correlations between the impurity potentials vS l > in different layers. The total scattering 
rates are defined by 



1 



w - - 



1 



(111) 



where the symbol (...) stands for the angular average. The disorder correlations between the layers are described by 
the characteristic rate 



1 



712 - T 



(112) 



where 1/r = [1/th + l/r22]/2. The time scale r g is a characteristic scale on which carriers in the two layers start 
"feeling" the difference between the impurity potentials and u^ 2 K The potentials in the two layers are strongly 
correlated when r g 3> r. One might expect that for realistic systems the situation of moderately correlated potentials, 
T g ~ r ~ Ti2, is typically realized. Weakly correlated potentials (712 3> t) yield r g <C r. Below we assume that 
disorder is sufficiently short-ranged and do not distinguish between the total and transport scattering rates for the 
estimates. 

We start from the ballistic regime Tt 3> 1. The correlated disorder affects the drag in a way similar to the third- 
order drag. With correlated disorder, one can include an interlayer disorder line w\ 2 into the inelastic scattering 

(3) 

amplitude. In the ballistic p D drag we had one amplitude M2 with two interaction lines and one with a single wave 
line (Mi). The corresponding drag rate contains 2Re[Mi(M2)*] oc a 3 . Now one can form the second-order scattering 
amplitude Mi using one interaction line (a) and one interlayer-disorder line, which introduces a factor (TT12) -1 . This 
gives 



1 



,-corr 
T D 



(113) 



and 



corr 

Pd 



Tt 12 



(114) 



which overcomes the third-order drag ~ a 3 for I/712 > aT. This happens in the perturbative regime (1/r > a 2 T, 
assuming correlated disorder, 712 ~ r), where the correlated-disorder contribution can be calculated diagrammatically. 
Similarly to a D , the corresponding diagram involves two four- leg vertices (hence finite drag at the Dirac point /i = 0), 
but now connected in all possible ways by two interaction lines and one disorder line Wi2- 

The general expression for the drag resistivity in the ballistic regime, including both third-order and correlated- 
disorder drag rates for equal layers has the form: 



h 1 



PD 



NT 



a 2 > 




Tt\i) 




H 2 N 


a 2 T + 


J^2 



a 4 T 2 - 



a 2 ^+a 3 T- 



c? y 
tti 2 j 



2lL_ 
J^2 



a 3 T 



Tt\i) 



(115) 



Exactly at the Dirac point it reduces to: 



Pd(h = 0) rC 



/ 1 



\Tti2 



(116) 



Let us now analyze the role of correlated disorder in the diffusive regime Tt <C 1. Again, we assume the absence of 
localization at the Dirac point (see Section 3). The drag resistivity for the case of correlated disorder was calculated 
in the diffusive regime in Ref. 



34 



It is dominated by the Maki-Thompson diagram with an interlayer Cooper 
propagator. It is worth noting that any difference in the disorder potentials (as well as in chemical potentials of the 
layers) leads to a finite gap in these propagators given by l/r g . The main result of Ref. 34) is as follows: 



corr 

Pd 



■ in 



Tt v t 9 



(117) 
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at Tg 1 «T« r _1 , and 



(Tt ) 2 

corr \ 91 /i i o\ 

P ° ~ 2 9^+He r g) r (118) 



at lower temperatures T <C t~ 1 . 

In graphene near the Dirac point, for small interlayer distance nd <C 1 the interlayer interaction constant in the 
Cooper channel is A12 ~ a. The Cooper channel cutoff energy is 60 = l/ r (the logarithm in the Cooper channel 
appears only for a constant density of states; in graphene in the diffusive regime this happens only for energies below 
1/t), the dimensionless conductance g ~ 1, and ~ 1/T. Substituting these values to Eqs. (jllTf) and (|118[) . we 
arrive at 

""" ~ (119) 

"°" ~ ? v-S$r)r T<<T - (120) 

These results are oc a 2 for realistic temperatures, Tt 3> exp(— 1/a). For a moderately correlated disorder r s ~ r, 
Eqs. (1TT4")) and ([T2T)1) then lead to 




which yields a maximum at T ~ 1/r in the temperature dependence of the drag resistivity at the charge neutrality 
point. For strongly correlated disorder potentials (r g 3> r), this maximum develops into a plateau between r" 1 <C 



B. Long-range correlations: correlated macroscopic inhomogeneities 

Let us now analyze within our kinetic-equation framework the model of correlated macroscopic spatial fluctuations 
Sfii in chemical potentials of the two layers (35j . characterized by the correlation function 

F^\r - r') = (8^(r)6^(r')) ? 0. (122) 

We restrict ourselves to the ballistic regime Tt ;§> 1. Assuming the spatial scale of the fluctuations to be much larger 
than all characteristic scales related to the particle scattering, vr ee , vtd, and i>t, we solve the hydrodynamic equations 
locally, yielding Eq. (pjTj) with local values of the chemical potentials encoded in functions ~ fi a ^/T, as well as 

in the local drag rate 

1 2 w MiOWr) 



r D {r) T 

On the other hand, since the coefficients Bq ~ 1 and B2 ~ 1, as well as the transport electron-electron scattering 
rate r" 1 ~ a 2 T are finite at the neutrality point, we can neglect the fluctuations of \ii in these quantities. Exactly 
at the Dirac point ^x,2 — ; assuming that the fluctuations of chemical potentials are weak (the precise condition 
is established below), we can further neglect the £?i-terms in the denominator of Eq. ([ST]) , yielding for the "local 
resistivity" 



Pd(v) 



2B ? t 



; 2 NBfiB^T 



1 Bl{r)B\{r) 



TT D {r) T? e T t 



b 

ee ' ee 



■ ^S^(r)5^(r) f^ + a ^N 2 ). ( \ 2.'! ) 



Averaging this expression over the small fluctuations of the correlated chemical potentials 35 1. we arrive at the 
correction to the universal third-order result, p%\n — 0) ~ (h/e 2 )a 3 , 

ApD ( M = 0) ~ » ^ I (1 + a 2 7VTr) ~ 2 " 2 1 (124) 

e 1 e 1 \a 2 a 2 iVT<-<T. 
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We see that in the Coulomb-dominated transport regime, this correction is dominated by the fluctuations in B\, 
whereas in the disorder-dominated (perturbative) regime, the main role is played by a locally finite drag rate. 
Finally, in the ultraclean limit 

- ^a 2 NFi 0) /T, (125) 

T 

we can neglect 1 jr in the denominator of the local drag resistivity given by Eq. (1611) , yielding a natural analog of 
Eq. CE3: 

Ap*(A. = 0)(r) ~ . (126) 

In particular, for perfectly correlated chemical potentials, <fyti(r) = J/i 2 (r), the fluctuations drops out from Eq. (|126[) 
and the local resistivity turns out to be independent of r. In a more general case, the averaging over fluctuations 
becomes nontrivial, but this can only affect the numerical prefactor in the final result. Thus, the correlated large-scale 
fluctuations of the chemical potentials in the layers in effect shift the curve 1 in Fig. [5] upwards, extending the validity 
of the fully equilibrated result, 




(127) 



to the case of finite disorder, Eq. (|125|) . at the Dirac point. This implies that in the case of correlated inhomogeneities 
the disorder-induced dip in the lower left panel of Fig. 1 develops only for sufficiently strong disorder. 
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